library(tidyverse)
library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)
library(readxl)
Greenhouse traits
Field traits
“H1: Pairing field results with a greenhouse 15 N experiment will allow us to relate species responses to their ability to uptake organic N”
“In addition, we will conduct a greenhouse experiment to isolate the mechanism (mineralized N vs organic N uptake rates ) by which compost may affect species composition”
Workflow plan:
overlay species trait data as vectors on the NMDS
only greenhouse traits available (only sla data from field)
use PC scores as composite trait axes OR
average the traits for each species
# # read in greenhouse trait data
# traits.gh <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=ya1ck6p7&dl=1")
#
# # read in field trait data
# # download.file("https://www.dropbox.com/scl/fi/trqrnz54dkjoa9do4tle5/leaf-area-data_NG.xlsx?rlkey=s0kic2a5yqknikc4icek8y7rn&st=dl4x634h&dl=1",destfile = "field_sla.xlsx")
# traits.sla.field <-
# read_excel("field_sla.xlsx",col_names = TRUE,col_types = "text")
#
# rda(traits.gh, scale = TRUE)
Trait data is identical. Will clean data in
# dropbox.main <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=0h1bs53j&dl=1")
#
# dropbox.gh <- read.csv("https://www.dropbox.com/scl/fi/67aaubdydiy4ixffmwvfe/GreenhouseTraits.csv?rlkey=rm2eqeqome7hgoedhnal0gx6d&st=vkym6evx&dl=1")
#
# all.equal(dropbox.main, dropbox.gh)
# identical(dropbox.main, dropbox.gh)
traits.cleaned <- read.csv("data/GreenhouseTraits_corrected_20240212.csv") %>%
mutate(
Fresh.leaf.mass..g. = as.numeric(ifelse(
Fresh.leaf.mass..g. == "Skipped accidentally", NA, Fresh.leaf.mass..g.
)),
Dry.leaf.mass..g. = as.numeric(ifelse(
Dry.leaf.mass..g. == "No sample", NA, Dry.leaf.mass..g.
)),
Root.dry.biomass..g. = as.numeric(ifelse(
Root.dry.biomass..g. == "roots lost", NA, Root.dry.biomass..g.
)),
Root.volume..cm3. = as.numeric(ifelse(
Root.volume..cm3. == "not scanned", NA, Root.volume..cm3.
))
)
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `Root.dry.biomass..g. = as.numeric(ifelse(Root.dry.biomass..g. == "roots lost", NA, Root.dry.biomass..g.))`.
Caused by warning:
! NAs introduced by coercion
# make trait species ID match matrix data (notes shown below) -------
traits.cleaned <- traits.cleaned %>%
mutate(
ID = ifelse(
ID == "Agoseris", "AGSP", ifelse(
ID == "AVEBAR", "AVBA", ifelse(
ID == "BRDR", "BRDI", ID
)
)
)
)
# summary(traits.cleaned)
traits.means <- traits.cleaned %>%
group_by(ID) %>%
summarise(
across(4:(ncol(traits.gh) - 1), list(mean = ~mean(.x, na.rm = TRUE)
))) %>%
rename(
Height = Height..cm._mean,
FreshLeafMass = Fresh.leaf.mass..g._mean,
DryLeafMass = Dry.leaf.mass..g._mean,
LDMC = LDMC_mean, # Leaf dry matter content (LDMC, the ratio of leaf dry mass to fresh mass)
LeafArea = Leaf.Area..cm2._mean,
SLA = SLA..cm2.g._mean,
ShootDryBiomass = Shoot.dry.biomass..g._mean,
RootDryBiomass = Root.dry.biomass..g._mean,
TotalBiomass = Total.biomass..g._mean,
RMF = RMF_mean, # Root mass fraction (RMF) is a plant trait that measures the proportion of a plant's dry mass that is in its roots
RootVol = Root.volume..cm3._mean,
RootDensity = Root.density..g.cm3._mean,
CoarseRootDiameter = Coarse.root.diameter..mm._mean,
RootLength = Length..mm._mean,
FineRootLength = Fine.root.length..mm._mean,
CoarseRootLength = Coarse.root.length..mm._mean,
CoarseRootSpecLength = Coarse.root.specific.length..cm.g._mean,
FineRootSpecLength = Fine.root.specific.length..cm.g._mean,
PropFineRoots = Proportion.fine.roots_mean
)
# make site by species matrix match traits species
matrix.names <- data.frame(colnames(matrix.bray)) #species comp (75 count)
trait.names <- data.frame(traits.means$ID) #trait species (80 count)
intersect.names <- intersect(colnames(matrix.bray), traits.means$ID) # same species (43 count)
setdiff(colnames(matrix.bray), traits.means$ID)
[1] "ASSP" "AST1" "AST2" "AST3" "CIQU" "CRTI" "CYEC" "DACA" "DIMU" "GAGR" "GAPA" "GAPH" "HYSP" "JUBU" "LYHY" "MAD1" "MAD2" "NAPU" "PAVI" "SABI1"
[21] "SABI2" "SIGA" "THGR" "TRSP" "UNBU" "UNF1" "UNF3" "UNF4" "UNF8" "UNFR" "UNGR" "XAST"
setdiff(traits.means$ID, colnames(matrix.bray))
[1] "ACMI" "ANARp" "AVFA" "BRCA" "BRHOp" "BRNIf" "BRNIp" "CLPUf" "CLPUp" "CYDA" "Crassula" "ESCA" "FEMI"
[14] "FEMY" "GITRf" "GITRp" "HYGL" "LACA" "LENIf" "LENIp" "LOPU" "LUBI" "Leontodon" "MAELf" "MAELp" "MICA"
[27] "NEMA" "PLERf" "PLERp" "RUPU" "Sanicula" "TRCI" "TRHIpi" "TRWIf" "TRWIp" "TRWIpi" "VIVA"
# check matrix codes to see if same species are named different things
unique.matrix.species <- cover.dat %>%
dplyr::select(code4, species) %>%
summarise(
code4 = unique(code4),
species = unique(species)
)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.
unique.trait.species <- traits.cleaned %>%
dplyr::select(Taxon, ID) %>%
distinct(ID, .keep_all = TRUE)
# write traits.means csv ----
# write.csv(traits.means, "data/traits/all_sp_trait_means.csv")
# # all years -------------
# # make site by species matrix match traits species
# matrix.names <- data.frame(colnames(matrix.bray)) #species comp (75 count)
# trait.names <- data.frame(traits.means$ID) #trait species (80 count)
#
# intersect.names <- intersect(colnames(matrix.bray), traits.means$ID) # same species (43 count)
# setdiff(colnames(matrix.bray), traits.means$ID)
# setdiff(traits.means$ID, colnames(matrix.bray))
#
# # check matrix codes to see if same species are named different things
# unique.matrix.species <- cover.dat %>%
# dplyr::select(code4, species) %>%
# summarise(
# code4 = unique(code4),
# species = unique(species)
# )
#
# unique.trait.species <- traits.cleaned %>%
# dplyr::select(Taxon, ID) %>%
# distinct(ID, .keep_all = TRUE)
#
# # 2019 -------------
# data.frame(colnames(matrix.bray.2019)) #species comp (75 count)
# data.frame(colnames(matrix.bray.2020)) #species comp (75 count)
# data.frame(colnames(matrix.bray.2021)) #species comp (75 count)
change Trait name: ‘Agoseris’ –> ‘Agoseris unknown’; assuming the same as ‘AGSP’ –> ‘Agoseris sp.’ in matrix data (change to matrix name)
What’s the difference between:
‘ANAR’ and ‘ANARp’ in trait data?
‘BRHO’ and ‘BRHOp’ in trait data?
‘BRNIf’ and ‘BRNIp’ in trait data? (not in matrix data)
‘CLPUf’ and ‘CLPUp’ in trait data? (not in matrix data)
‘GITRf’ and ‘GITRp’ in trait data? (not in matrix data)
‘LENIf’ and ‘LENIp’ in trait data? (not in matrix data)
‘TRHI’ and ‘TRHIpi’ (same sp name) in trait data? (TRHI in matrix data)
‘TRWIf’ and ‘TRWIp’ and ‘TRWIpi’ in trait data? (not in matrix data)
‘Crassula’ and ‘Crassula tillaea’ (could they be the same?)
‘Sanicula’ (unknown) and ‘SABI1’ or ‘SABI2’ (could they be the same?)
‘MAELf’ (trait data) and ‘MAD1’ (Madia sp. 1 (“tall tarweed”)’ (could they be the same?)
‘MAELp’ (trait data) and ‘MAD2’ (Madia sp. 1 (“small tarweed”)’ (could they be the same?)
‘Trifolium eriocephalum’ (TRER) and ‘Triphysaria eriantha’ (TRER) (could they be the same?)
change ‘AVEBAR’ (trait data) to ‘AVBA’ (matrix data)
change BRDR –> BRDI
Aboveground Traits
Height
Leaf
FreshLeafMass
DryLeafMass
LDMC
LeafArea
SLA
ShootDryBiomass
Belowground Traits
RootDryBiomass
RMF
RootVol
RootDensity
RootLength
PropFineRoots
Coarse Roots
CoarseRootDiameter
CoarseRootLength
CoarseRootSpecLength
Fine Roots
FineRootLength
FineRootSpecLength
TotalBiomass
Stress:
0.2011127
# subset nmds to match species names
nmds.comp.trait <- metaMDS(matrix.bray %>%
dplyr::select(intersect.names),
k=2, trymax = 25)
Run 0 stress 0.2060005
Run 1 stress 0.2019906
... New best solution
... Procrustes: rmse 0.07567441 max resid 0.2449078
Run 2 stress 0.2007527
... New best solution
... Procrustes: rmse 0.08985478 max resid 0.2109131
Run 3 stress 0.2094711
Run 4 stress 0.2107847
Run 5 stress 0.2089838
Run 6 stress 0.2078775
Run 7 stress 0.2027916
Run 8 stress 0.2173911
Run 9 stress 0.2059157
Run 10 stress 0.2220001
Run 11 stress 0.2079152
Run 12 stress 0.2124366
Run 13 stress 0.2050097
Run 14 stress 0.202904
Run 15 stress 0.2027187
Run 16 stress 0.2032121
Run 17 stress 0.2179617
Run 18 stress 0.2028023
Run 19 stress 0.2012207
... Procrustes: rmse 0.02040594 max resid 0.1298012
Run 20 stress 0.2051096
Run 21 stress 0.2080723
Run 22 stress 0.201363
Run 23 stress 0.2029037
Run 24 stress 0.2173911
Run 25 stress 0.2079715
*** Best solution was not repeated -- monoMDS stopping criteria:
22: stress ratio > sratmax
3: scale factor of the gradient < sfgrmin
nmds.comp.trait # pretty not great stress
Call:
metaMDS(comm = matrix.bray %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.2007527
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 2 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray %>% dplyr::select(intersect.names)’
stressplot(nmds.comp.trait)
plot(nmds.comp.trait, type="t")
nmds1 <- as.data.frame(scores(nmds.comp.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.comp.trait, choices=c(2), display=c("sites")))
nmds_dat_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# subset matrix and trait data to match names
traits.means.comp <- traits.means %>%
filter(ID %in% intersect.names)
# traits.means.comp2 <- traits.means.comp %>%
# dplyr::select(-ID)
rownames(traits.means.comp2) <- traits.means.comp$ID
Warning: Setting row names on a tibble is deprecated.
# trait vectors onto the NMDS
# envfit(nmds.comp.trait, traits.means.comp, permutations = 999)
# example from stack overflow: envfit(ord ~ var1 + var2 + var3, dune.spec, display="sp"), 'ord' is the nmds model, 'var1' is the trait1 i think?
colnames(traits.means.comp)
[1] "ID" "Height" "FreshLeafMass" "DryLeafMass" "LDMC" "LeafArea" "SLA"
[8] "ShootDryBiomass" "RootDryBiomass" "TotalBiomass" "RMF" "RootVol" "RootDensity" "CoarseRootDiameter"
[15] "RootLength" "FineRootLength" "CoarseRootLength" "CoarseRootSpecLength" "FineRootSpecLength" "PropFineRoots"
vectors.nmds.comp.trait <- envfit(nmds.comp.trait ~ Height +
LeafArea +
RMF, traits.means.comp, display = "sp") # ONLY Height, Leaf area, and RMF significant
#TEMPLATE:
# envfit(nmds.block.trait ~ Height +
# FreshLeafMass +
# DryLeafMass +
# LDMC +
# LeafArea +
# SLA +
# ShootDryBiomass +
# RootDryBiomass +
# TotalBiomass +
# RMF +
# RootVol +
# RootDensity +
# CoarseRootDiameter +
# RootLength +
# FineRootLength +
# CoarseRootLength +
# CoarseRootSpecLength +
# FineRootSpecLength +
# PropFineRoots, traits.means.comp, display = "sp", na.rm = TRUE)
# is it different when you only do certain traits (above- vs. below-ground)? or can you just throw everything in there and it will tell you what is significant and you can subset from there?
# Looks like different subsets of traits put in are the same, so I shall proceed with the mega models and verify with Lauren
trait_vectors <- as.data.frame(vectors.nmds.comp.trait$vectors$arrows)
NMDS model: only Height, Leaf area, and RMF () found significant.
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year")+
geom_segment(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)+
geom_segment(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Grazing History") +
geom_segment(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
# ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = as.factor(grazing_hist), fill = as.factor(yr))) +
# geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
# stat_ellipse(aes(color = as.factor(yr)),
# geom = "polygon",
# fill = NA,
# size = 1) +
# labs(title = "Grouped by Grazing History & Year") +
# scale_color_manual(
# values = c(
# "low" = "cyan",
# "high" = "brown",
# "2019" = "salmon",
# "2020" = "limegreen",
# "2021" = "cornflowerblue")
# ) +
# geom_segment(
# data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
# aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
# inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
# arrow = arrow(length = unit(0.2, "cm")),
# color = "red",
# size = 1
# ) +
# geom_text(
# data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
# aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)),
# inherit.aes = FALSE, # Avoids unwanted inheritance
# hjust = -0.2,
# vjust = -0.2,
# color = "red",
# size = 5
# )
#
# ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
# geom_point(shape = 21, size = 4, stroke = 1) +
# stat_ellipse(aes(color = as.factor(nut_trt)),
# geom = "polygon",
# fill = NA, # Makes the ellipse transparent inside
# size = 1) +
# labs(title = "Grouped by Nutrient Treatment & Year") +
# scale_color_manual(
# values = c(
# "control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
# "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
# "+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
# "wet" = rgb(0, 0, 255, maxColorValue = 255),
# "dry" = rgb(183, 65, 14, maxColorValue = 255)
# )
# ) +
# geom_segment(
# data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
# aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
# inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
# arrow = arrow(length = unit(0.2, "cm")),
# color = "red",
# size = 1
# ) +
# geom_text(
# data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
# aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)),
# inherit.aes = FALSE, # Avoids unwanted inheritance
# hjust = -0.2,
# vjust = -0.2,
# color = "red",
# size = 5
# )
Stress:
0.2353937
# subset nmds to match species names
nmds.block.trait <- metaMDS(matrix.block.bray %>%
dplyr::select(intersect.names), , k=2, trymax = 25)
Run 0 stress 0.2353937
Run 1 stress 0.250124
Run 2 stress 0.2353908
... New best solution
... Procrustes: rmse 0.002559202 max resid 0.02324147
Run 3 stress 0.2386994
Run 4 stress 0.2487719
Run 5 stress 0.2453858
Run 6 stress 0.2574031
Run 7 stress 0.2354028
... Procrustes: rmse 0.002945858 max resid 0.02371427
Run 8 stress 0.2563643
Run 9 stress 0.2460121
Run 10 stress 0.2644327
Run 11 stress 0.2560224
Run 12 stress 0.246987
Run 13 stress 0.2354028
... Procrustes: rmse 0.002938524 max resid 0.02369694
Run 14 stress 0.2423148
Run 15 stress 0.2463114
Run 16 stress 0.2354028
... Procrustes: rmse 0.002925757 max resid 0.02366134
Run 17 stress 0.2353936
... Procrustes: rmse 0.002546433 max resid 0.02319281
Run 18 stress 0.2526558
Run 19 stress 0.2499386
Run 20 stress 0.2480986
Run 21 stress 0.2497376
Run 22 stress 0.2409873
Run 23 stress 0.2385602
Run 24 stress 0.2521044
Run 25 stress 0.2378706
*** Best solution was not repeated -- monoMDS stopping criteria:
24: stress ratio > sratmax
1: scale factor of the gradient < sfgrmin
nmds.block.trait # pretty not great stress
Call:
metaMDS(comm = matrix.block.bray %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.block.bray %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.2353908
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 2 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.block.bray %>% dplyr::select(intersect.names)’
nmds1 <- as.data.frame(scores(nmds.block.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.block.trait, choices=c(2), display=c("sites")))
nmds_dat_block_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# trait data
traits.means.comp
# trait vectors onto the NMDS
vectors.nmds.block.trait <-
envfit(nmds.block.trait ~ Height +
DryLeafMass +
LeafArea +
RootDryBiomass +
RMF +
RootDensity, traits.means.comp, display = "sp")
# ONLY Height, DryLeafMass, LeafArea, RMF significant (RootDryBiomass, RootDensity marginal)
trait_vectors <- as.data.frame(vectors.nmds.block.trait$vectors$arrows)
ggplot(nmds_dat_block_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment - Blocks w Traits") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_block_trait, aes(NMDS1, NMDS2, fill = ppt_trt, color = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment - Blocks w Traits") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
NA
NA
NA
NA
Stress:
0.1879671
# subset nmds to match species names
nmds.low.trait <- metaMDS(matrix.bray.low %>%
dplyr::select(intersect.names), , k=2, trymax = 25)
Run 0 stress 0.2032274
Run 1 stress 0.1911123
... New best solution
... Procrustes: rmse 0.07305776 max resid 0.1904385
Run 2 stress 0.2933553
Run 3 stress 0.1890672
... New best solution
... Procrustes: rmse 0.1307018 max resid 0.279877
Run 4 stress 0.2104987
Run 5 stress 0.208207
Run 6 stress 0.2170339
Run 7 stress 0.2019759
Run 8 stress 0.1911123
Run 9 stress 0.2124695
Run 10 stress 0.2032282
Run 11 stress 0.2020341
Run 12 stress 0.1890672
... Procrustes: rmse 3.79248e-06 max resid 8.350358e-06
... Similar to previous best
Run 13 stress 0.2138031
Run 14 stress 0.3881609
Run 15 stress 0.2019912
Run 16 stress 0.1915469
Run 17 stress 0.2112419
Run 18 stress 0.2084351
Run 19 stress 0.1911123
Run 20 stress 0.2105596
*** Best solution repeated 1 times
nmds.low.trait
Call:
metaMDS(comm = matrix.bray.low %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.low %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.1890672
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 3 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.low %>% dplyr::select(intersect.names)’
nmds1 <- as.data.frame(scores(nmds.low.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.low.trait, choices=c(2), display=c("sites")))
nmds_dat_low_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# trait data: traits.means.comp
# trait vectors onto the NMDS
vectors.nmds.low.trait <-
envfit(nmds.low.trait ~ FreshLeafMass +
RMF, traits.means.comp, display = "sp", na.rm = TRUE)
# ONLY FreshLeafMass (RMF marginal)
trait_vectors <- as.data.frame(vectors.nmds.low.trait$vectors$arrows)
ggplot(nmds_dat_low_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
labs(title = "2019-2021 (LOW) w Traits - Grouping by Nutrient Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_low_trait, aes(NMDS1, NMDS2, fill = ppt_trt, color = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
labs(title = "2019-2021 (LOW) w Traits - Grouping by Precipitation Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
Stress:
0.06290431
# subset nmds to match species names
nmds.2019.trait <- metaMDS(matrix.bray.2019 %>%
dplyr::select(intersect.names), , k=2, trymax = 25)
Run 0 stress 0.06471496
Run 1 stress 0.06931106
Run 2 stress 0.06874971
Run 3 stress 0.3616162
Run 4 stress 0.06264873
... New best solution
... Procrustes: rmse 0.06539037 max resid 0.2367651
Run 5 stress 0.09229892
Run 6 stress 0.06471497
Run 7 stress 0.06290424
... Procrustes: rmse 0.005564852 max resid 0.0187599
Run 8 stress 0.06816428
Run 9 stress 0.06901458
Run 10 stress 0.06471496
Run 11 stress 0.07109362
Run 12 stress 0.06459268
Run 13 stress 0.08080796
Run 14 stress 0.07203292
Run 15 stress 0.06459267
Run 16 stress 0.07780829
Run 17 stress 0.07518585
Run 18 stress 0.06264873
... Procrustes: rmse 1.714225e-05 max resid 4.819798e-05
... Similar to previous best
Run 19 stress 0.06931106
Run 20 stress 0.07273526
*** Best solution repeated 1 times
nmds.2019.trait
Call:
metaMDS(comm = matrix.bray.2019 %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2019 %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.06264873
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 4 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2019 %>% dplyr::select(intersect.names)’
nmds1 <- as.data.frame(scores(nmds.2019.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2019.trait, choices=c(2), display=c("sites")))
nmds_dat_2019_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# trait data: traits.means.comp
# trait vectors onto the NMDS
vectors.nmds.2019.trait <-
envfit(nmds.2019.trait ~ Height +
DryLeafMass +
LeafArea +
RMF +
RootDensity +
CoarseRootDiameter, traits.means.comp, display = "sp", na.rm = TRUE)
ggplot(nmds_dat_2019_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2019 w Traits - Grouping by Nutrient Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_2019_trait, aes(NMDS1, NMDS2, fill = ppt_trt, color = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2019 w Traits - Grouping by Precipitation Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
Stress:
0.1512864
# subset nmds to match species names
nmds.2020.trait <- metaMDS(matrix.bray.2020 %>%
dplyr::select(intersect.names), , k=2, trymax = 25)
Run 0 stress 0.1512864
Run 1 stress 0.16382
Run 2 stress 0.1740134
Run 3 stress 0.1512864
... New best solution
... Procrustes: rmse 1.160478e-05 max resid 2.59613e-05
... Similar to previous best
Run 4 stress 0.151653
... Procrustes: rmse 0.01975231 max resid 0.06777473
Run 5 stress 0.16382
Run 6 stress 0.1601231
Run 7 stress 0.151336
... Procrustes: rmse 0.01055512 max resid 0.03081296
Run 8 stress 0.1600829
Run 9 stress 0.1922664
Run 10 stress 0.171301
Run 11 stress 0.1512864
... New best solution
... Procrustes: rmse 1.626593e-05 max resid 3.63373e-05
... Similar to previous best
Run 12 stress 0.151336
... Procrustes: rmse 0.0105455 max resid 0.03079015
Run 13 stress 0.1600829
Run 14 stress 0.163364
Run 15 stress 0.1512864
... Procrustes: rmse 8.079558e-06 max resid 2.224799e-05
... Similar to previous best
Run 16 stress 0.1740133
Run 17 stress 0.1516619
... Procrustes: rmse 0.01548936 max resid 0.05232438
Run 18 stress 0.1714389
Run 19 stress 0.175096
Run 20 stress 0.16382
*** Best solution repeated 2 times
nmds.2020.trait
Call:
metaMDS(comm = matrix.bray.2020 %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2020 %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.1512864
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 11 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2020 %>% dplyr::select(intersect.names)’
nmds1 <- as.data.frame(scores(nmds.2020.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2020.trait, choices=c(2), display=c("sites")))
nmds_dat_2020_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# trait data: traits.means.comp
# trait vectors onto the NMDS
vectors.nmds.2020.trait <-
envfit(nmds.2020.trait ~ FreshLeafMass +
LeafArea, traits.means.comp, display = "sp", na.rm = TRUE)
# ONLY FreshLeafMass, LeafArea
trait_vectors <- as.data.frame(vectors.nmds.2020.trait$vectors$arrows)
ggplot(nmds_dat_2020_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 w Traits - Grouping by Nutrient Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_2020_trait, aes(NMDS1, NMDS2, fill = ppt_trt, color = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 w Traits - Grouping by Precipitation Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
Stress:
0.123061
# subset nmds to match species names
nmds.2021.trait <- metaMDS(matrix.bray.2021 %>%
dplyr::select(intersect.names), , k=2, trymax = 25)
Run 0 stress 0.1234883
Run 1 stress 0.123061
... New best solution
... Procrustes: rmse 0.04524088 max resid 0.133887
Run 2 stress 0.1556472
Run 3 stress 0.1556472
Run 4 stress 0.1510891
Run 5 stress 0.133126
Run 6 stress 0.1412774
Run 7 stress 0.1234883
... Procrustes: rmse 0.04524054 max resid 0.133303
Run 8 stress 0.1234883
... Procrustes: rmse 0.04524098 max resid 0.13329
Run 9 stress 0.1331269
Run 10 stress 0.1412774
Run 11 stress 0.1331264
Run 12 stress 0.1329969
Run 13 stress 0.123061
... Procrustes: rmse 7.588273e-06 max resid 2.28136e-05
... Similar to previous best
Run 14 stress 0.1412774
Run 15 stress 0.1412774
Run 16 stress 0.1234883
... Procrustes: rmse 0.04524105 max resid 0.1332899
Run 17 stress 0.1234883
... Procrustes: rmse 0.04524004 max resid 0.1333063
Run 18 stress 0.1790971
Run 19 stress 0.1513012
Run 20 stress 0.123061
... Procrustes: rmse 2.514183e-06 max resid 7.876099e-06
... Similar to previous best
*** Best solution repeated 2 times
nmds.2021.trait
Call:
metaMDS(comm = matrix.bray.2021 %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2021 %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.123061
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 1 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2021 %>% dplyr::select(intersect.names)’
nmds1 <- as.data.frame(scores(nmds.2021.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2021.trait, choices=c(2), display=c("sites")))
nmds_dat_2021_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# trait data: traits.means.comp
rownames(traits.means.comp) <- traits.means.comp$ID
Warning: Setting row names on a tibble is deprecated.
# trait vectors onto the NMDS
vectors.nmds.2021.trait <-
envfit(nmds.2021.trait ~ Height +
FreshLeafMass +
LeafArea, na.omit(traits.means.comp), display = "sp", na.rm = TRUE)
# (Height, FreshLeafMass, LeafArea marginal)
trait_vectors <- as.data.frame(vectors.nmds.2021.trait$vectors$arrows)
ggplot(nmds_dat_2021_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 w Traits - Grouping by Nutrient Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_2021_trait, aes(NMDS1, NMDS2, fill = ppt_trt, color = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 w Traits - Grouping by Precipitation Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
Stress:
0.170799
# subset nmds to match species names
nmds.2021.block.trait <- metaMDS(matrix.bray.block.2021 %>%
dplyr::select(intersect.names), , k=2, trymax = 25)
Run 0 stress 0.170799
Run 1 stress 0.2052476
Run 2 stress 0.170799
... New best solution
... Procrustes: rmse 9.313132e-06 max resid 3.473702e-05
... Similar to previous best
Run 3 stress 0.2060774
Run 4 stress 0.1873003
Run 5 stress 0.2046024
Run 6 stress 0.206467
Run 7 stress 0.2180152
Run 8 stress 0.1709819
... Procrustes: rmse 0.007084111 max resid 0.03059497
Run 9 stress 0.1797966
Run 10 stress 0.170799
... Procrustes: rmse 8.637585e-06 max resid 3.219545e-05
... Similar to previous best
Run 11 stress 0.1797966
Run 12 stress 0.1914994
Run 13 stress 0.1951365
Run 14 stress 0.1996107
Run 15 stress 0.196785
Run 16 stress 0.170799
... New best solution
... Procrustes: rmse 6.097064e-06 max resid 2.279654e-05
... Similar to previous best
Run 17 stress 0.214524
Run 18 stress 0.1709819
... Procrustes: rmse 0.007078194 max resid 0.03057796
Run 19 stress 0.189269
Run 20 stress 0.1939681
*** Best solution repeated 1 times
nmds.2021.block.trait
Call:
metaMDS(comm = matrix.bray.block.2021 %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.block.2021 %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.170799
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 16 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.block.2021 %>% dplyr::select(intersect.names)’
nmds1 <- as.data.frame(scores(nmds.2021.block.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2021.block.trait, choices=c(2), display=c("sites")))
nmds_dat_2021_block_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# trait data: traits.means.comp
# trait vectors onto the NMDS
vectors.nmds.2021.block.trait <-
envfit(nmds.2021.block.trait ~ DryLeafMass +
LeafArea, traits.means.comp, display = "sp", na.rm = TRUE)
# ONLY DryLeafMass, LeafArea
trait_vectors <- as.data.frame(vectors.nmds.2021.block.trait$vectors$arrows)
ggplot(nmds_dat_2021_block_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block) w Traits - Grouping by Nutrient Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_2021_block_trait, aes(NMDS1, NMDS2, fill = ppt_trt, color = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block) w Traits - Grouping by Precipitation Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
Stress:
0.1299711
# subset nmds to match species names
nmds.2021.block.low.trait <- metaMDS(matrix.bray.block.2021.low %>%
dplyr::select(intersect.names), , k=2, trymax = 25)
Run 0 stress 0.1299711
Run 1 stress 0.1299712
... Procrustes: rmse 5.590124e-05 max resid 0.0001635787
... Similar to previous best
Run 2 stress 0.1299712
... Procrustes: rmse 4.398462e-05 max resid 0.0001319322
... Similar to previous best
Run 3 stress 0.1627521
Run 4 stress 0.1299711
... Procrustes: rmse 1.009167e-05 max resid 3.112151e-05
... Similar to previous best
Run 5 stress 0.162752
Run 6 stress 0.1299711
... Procrustes: rmse 1.406545e-05 max resid 4.211068e-05
... Similar to previous best
Run 7 stress 0.1299711
... Procrustes: rmse 7.661138e-06 max resid 1.674758e-05
... Similar to previous best
Run 8 stress 0.18467
Run 9 stress 0.2001228
Run 10 stress 0.1602444
Run 11 stress 0.1299712
... Procrustes: rmse 4.253368e-05 max resid 0.0001252648
... Similar to previous best
Run 12 stress 0.1431163
Run 13 stress 0.2111627
Run 14 stress 0.1431163
Run 15 stress 0.2104642
Run 16 stress 0.1299712
... Procrustes: rmse 3.767893e-05 max resid 0.0001028459
... Similar to previous best
Run 17 stress 0.1431163
Run 18 stress 0.1785202
Run 19 stress 0.1299711
... New best solution
... Procrustes: rmse 8.487148e-06 max resid 2.52772e-05
... Similar to previous best
Run 20 stress 0.1602445
*** Best solution repeated 1 times
nmds.2021.block.low.trait
Call:
metaMDS(comm = matrix.bray.block.2021.low %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.block.2021.low %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.1299711
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 19 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.block.2021.low %>% dplyr::select(intersect.names)’
nmds1 <- as.data.frame(scores(nmds.2021.block.low.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2021.block.low.trait, choices=c(2), display=c("sites")))
nmds_dat_2021_block_low_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# trait data: traits.means.comp
# trait vectors onto the NMDS
vectors.nmds.2021.block.low.trait <-
envfit(nmds.2021.block.low.trait ~ FreshLeafMass, traits.means.comp, display = "sp", na.rm = TRUE)
# ONLY FreshLeafMass
trait_vectors <- as.data.frame(vectors.nmds.2021.block.low.trait$vectors$arrows)
ggplot(nmds_dat_2021_block_low_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; LOW) w Trait - Grouping by Nutrient Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_2021_block_low_trait, aes(NMDS1, NMDS2, fill = ppt_trt, color = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; LOW) w Trait - Grouping by Precipitation Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = trait_vectors,
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = trait_vectors,
aes(x = NMDS1, y = NMDS2, label = rownames(trait_vectors)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
# subset nmds to match species names
nmds.2021.block.high.trait <- metaMDS(matrix.bray.block.2021.high %>%
dplyr::select(intersect.names), , k=2, trymax = 25)
Run 0 stress 0.1416462
Run 1 stress 0.1595099
Run 2 stress 0.1867893
Run 3 stress 0.2747088
Run 4 stress 0.1867893
Run 5 stress 0.1416462
... New best solution
... Procrustes: rmse 3.289538e-05 max resid 8.812783e-05
... Similar to previous best
Run 6 stress 0.1867893
Run 7 stress 0.1416461
... New best solution
... Procrustes: rmse 0.0004825076 max resid 0.00160539
... Similar to previous best
Run 8 stress 0.1416461
... Procrustes: rmse 0.0004023531 max resid 0.00133908
... Similar to previous best
Run 9 stress 0.1595099
Run 10 stress 0.1867893
Run 11 stress 0.15951
Run 12 stress 0.1867893
Run 13 stress 0.1416461
... Procrustes: rmse 7.442382e-05 max resid 0.000245897
... Similar to previous best
Run 14 stress 0.1416461
... Procrustes: rmse 4.391587e-05 max resid 0.000145225
... Similar to previous best
Run 15 stress 0.1595099
Run 16 stress 0.141646
... New best solution
... Procrustes: rmse 0.0002924508 max resid 0.0009736183
... Similar to previous best
Run 17 stress 0.1595099
Run 18 stress 0.1416461
... Procrustes: rmse 0.000208078 max resid 0.0006900715
... Similar to previous best
Run 19 stress 0.141646
... Procrustes: rmse 0.0002689142 max resid 0.0008954326
... Similar to previous best
Run 20 stress 0.1416463
... Procrustes: rmse 0.0002338219 max resid 0.0007719694
... Similar to previous best
*** Best solution repeated 4 times
nmds.2021.block.high.trait
Call:
metaMDS(comm = matrix.bray.block.2021.high %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.block.2021.high %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.141646
Stress type 1, weak ties
Best solution was repeated 4 times in 20 tries
The best solution was from try 16 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.block.2021.high %>% dplyr::select(intersect.names)’
nmds1 <- as.data.frame(scores(nmds.2021.block.high.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.2021.block.high.trait, choices=c(2), display=c("sites")))
nmds_dat_2021_block_high_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# trait data: traits.means.comp
# trait vectors onto the NMDS
vectors.nmds.2021.block.high.trait <-
envfit(nmds.2021.block.high.trait ~ Height +
FreshLeafMass +
DryLeafMass +
LDMC +
LeafArea +
SLA +
ShootDryBiomass +
RootDryBiomass +
TotalBiomass +
RMF +
RootVol +
RootDensity +
CoarseRootDiameter +
RootLength +
FineRootLength +
CoarseRootLength +
CoarseRootSpecLength +
FineRootSpecLength +
PropFineRoots, traits.means.comp, display = "sp", na.rm = TRUE)
# None
trait_vectors <- as.data.frame(vectors.nmds.2021.block.high.trait$vectors$arrows)
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; HIGH) - Grouping by Nutrient Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, fill = ppt_trt, color = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; HIGH) - Grouping by Precipitation Treatment") +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)